{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "D2J3nB-ZrRv1"
      },
      "source": [
        "##### Copyright 2018 The TensorFlow Authors.\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "9qDhTJmprPnm"
      },
      "outputs": [],
      "source": [
        "#@title Licensed under the Apache License, Version 2.0 (the \"License\"); { display-mode: \"form\" }\n",
        "# you may not use this file except in compliance with the License.\n",
        "# You may obtain a copy of the License at\n",
        "#\n",
        "# https://www.apache.org/licenses/LICENSE-2.0\n",
        "#\n",
        "# Unless required by applicable law or agreed to in writing, software\n",
        "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "# See the License for the specific language governing permissions and\n",
        "# limitations under the License."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "pfPtIQ3DdZ8r"
      },
      "source": [
        "# Generalized Linear Models\n",
        "\u003ctable class=\"tfo-notebook-buttons\" align=\"left\"\u003e\n",
        "  \u003ctd\u003e\n",
        "    \u003ca target=\"_blank\" href=\"https://colab.research.google.com/github/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Generalized_Linear_Models.ipynb\"\u003e\u003cimg src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" /\u003eRun in Google Colab\u003c/a\u003e\n",
        "  \u003c/td\u003e\n",
        "  \u003ctd\u003e\n",
        "    \u003ca target=\"_blank\" href=\"https://github.com/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Generalized_Linear_Models.ipynb\"\u003e\u003cimg src=\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\" /\u003eView source on GitHub\u003c/a\u003e\n",
        "  \u003c/td\u003e\n",
        "\u003c/table\u003e"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "EOfH1_F9YsOG"
      },
      "source": [
        "In this notebook we introduce Generalized Linear Models via a worked example.  We solve this example in two different ways using two algorithms for efficiently fitting GLMs in TensorFlow Probability: Fisher scoring for dense data, and coordinatewise proximal gradient descent for sparse data.  We compare the fitted coefficients to the true coefficients and, in the case of coordinatewise proximal gradient descent, to the output of R's similar `glmnet` algorithm.  Finally, we provide further mathematical details and derivations of several key properties of GLMs."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "rjsfQ6vLb5I0"
      },
      "source": [
        "# Background"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "TdMX-QKagFnY"
      },
      "source": [
        "A generalized linear model (GLM) is a linear model ($\\eta = x^\\top \\beta$) wrapped in a transformation (link function) and equipped with a response distribution from an exponential family.  The choice of link function and response distribution is very flexible, which lends great expressivity to GLMs.  The full details, including a sequential presentation of all the definitions and results building up to GLMs in unambiguous notation, are found in \"Derivation of GLM Facts\" below.  We summarize:\n",
        "\n",
        "In a GLM, a predictive distribution for the response variable $Y$ is associated with a vector of observed predictors $x$.  The distribution has the form:\n",
        "\n",
        "\\begin{align*}\n",
        "  p(y \\, |\\, x)\n",
        "\u0026=\n",
        "  m(y, \\phi) \\exp\\left(\\frac{\\theta\\, T(y) - A(\\theta)}{\\phi}\\right)\n",
        "\\\\\n",
        "  \\theta\n",
        "\u0026:=\n",
        "  h(\\eta)\n",
        "\\\\\n",
        "  \\eta\n",
        "\u0026:=\n",
        "  x^\\top \\beta\n",
        "\\end{align*}\n",
        "\n",
        "Here $\\beta$ are the parameters (\"weights\"), $\\phi$ a hyperparameter representing dispersion (\"variance\"), and $m$, $h$, $T$, $A$ are characterized by the user-specified model family.\n",
        "\n",
        "The mean of $Y$ depends on $x$ by composition of **linear response** $\\eta$ and (inverse) link function, i.e.:\n",
        "\n",
        "$$\n",
        "\\mu := g^{-1}(\\eta)\n",
        "$$\n",
        "\n",
        "where $g$ is the so-called **link function**.  In TFP  the choice of link function and model family are jointly specifed by a `tfp.glm.ExponentialFamily` subclass. Examples include:\n",
        "- `tfp.glm.Normal`, aka \"linear regression\"\n",
        "- `tfp.glm.Bernoulli`, aka \"logistic regression\"\n",
        "- `tfp.glm.Poisson`, aka \"Poisson regression\"\n",
        "- `tfp.glm.BernoulliNormalCDF`, aka \"probit regression\".\n",
        "\n",
        "TFP prefers to name model families according to the distribution over `Y` rather than the link function since `tfp.Distribution`s are already first-class citizens. If the `tfp.glm.ExponentialFamily` subclass name contains a second word, this indicates a [non-canonical link function](https://en.wikipedia.org/wiki/Generalized_linear_model#Link_function)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "1oGScpRnqH_b"
      },
      "source": [
        "GLMs have several remarkable properties which permit efficient implementation of the maximum likelihood estimator.  Chief among these properties are simple formulas for the gradient of the log-likelihood $\\ell$, and for the Fisher information matrix, which is the expected value of the Hessian of the negative log-likelihood under a re-sampling of the response under the same predictors.  I.e.:\n",
        "\n",
        "\\begin{align*}\n",
        "  \\nabla_\\beta\\, \\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{y})\n",
        "\u0026=\n",
        "  \\mathbf{x}^\\top\n",
        "  \\,\\text{diag}\\left(\\frac{\n",
        "      {\\textbf{Mean}_T}'(\\mathbf{x} \\beta)\n",
        "    }{\n",
        "      {\\textbf{Var}_T}(\\mathbf{x} \\beta)\n",
        "    }\\right)\n",
        "  \\left(\\mathbf{T}(\\mathbf{y}) - {\\textbf{Mean}_T}(\\mathbf{x} \\beta)\\right)\n",
        "\\\\\n",
        "  \\mathbb{E}_{Y_i \\sim \\text{GLM} | x_i} \\left[\n",
        "    \\nabla_\\beta^2\\, \\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{Y})\n",
        "  \\right]\n",
        "\u0026=\n",
        "  -\\mathbf{x}^\\top\n",
        "  \\,\\text{diag}\\left(\n",
        "    \\frac{\n",
        "      \\phi\\, {\\textbf{Mean}_T}'(\\mathbf{x} \\beta)^2\n",
        "    }{\n",
        "      {\\textbf{Var}_T}(\\mathbf{x} \\beta)\n",
        "    }\\right)\\,\n",
        "  \\mathbf{x}\n",
        "\\end{align*}\n",
        "\n",
        "where $\\mathbf{x}$ is the matrix whose $i$th row is the predictor vector for the $i$th data sample, and $\\mathbf{y}$ is vector whose $i$th coordinate is the observed response for the $i$th data sample.  Here (loosely speaking), ${\\text{Mean}_T}(\\eta) := \\mathbb{E}[T(Y)\\,|\\,\\eta]$ and ${\\text{Var}_T}(\\eta) := \\text{Var}[T(Y)\\,|\\,\\eta]$, and boldface denotes vectorization of these functions.  Full details of what distributions these expectations and variances are over can be found in \"Derivation of GLM Facts\" below."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "XuNDwfwBObKl"
      },
      "source": [
        "# An Example\n",
        "\n",
        "In this section we briefly describe and showcase two built-in GLM fitting algorithms in TensorFlow Probability: Fisher scoring (`tfp.glm.fit`) and coordinatewise proximal gradient descent (`tfp.glm.fit_sparse`)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "4phryMfsP4Sn"
      },
      "source": [
        "## Synthetic Data Set\n",
        "\n",
        "Let's pretend to load some training data set."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "DA2Rf9PPgMAD"
      },
      "outputs": [],
      "source": [
        "from __future__ import print_function\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import scipy\n",
        "import tensorflow as tf\n",
        "import tensorflow_probability as tfp\n",
        "tfd = tfp.distributions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "KEVnTz2hh9RN"
      },
      "outputs": [],
      "source": [
        "  def make_dataset(n, d, link, scale=1., dtype=np.float32):\n",
        "    model_coefficients = tfd.Uniform(\n",
        "        low=np.array(-1, dtype), high=np.array(1, dtype)).sample(\n",
        "            d, seed=42)\n",
        "    radius = np.sqrt(2.)\n",
        "    model_coefficients *= radius / tf.linalg.norm(model_coefficients)\n",
        "    mask = tf.random_shuffle(tf.range(d)) \u003c tf.to_int32(0.5 * tf.to_float(d))\n",
        "    model_coefficients = tf.where(mask, model_coefficients,\n",
        "                                  tf.zeros_like(model_coefficients))\n",
        "    model_matrix = tfd.Normal(\n",
        "        loc=np.array(0, dtype), scale=np.array(1, dtype)).sample(\n",
        "            [n, d], seed=43)\n",
        "    scale = tf.convert_to_tensor(scale, dtype)\n",
        "    linear_response = tf.matmul(model_matrix,\n",
        "                                model_coefficients[..., tf.newaxis])[..., 0]\n",
        "    if link == 'linear':\n",
        "      response = tfd.Normal(loc=linear_response, scale=scale).sample(seed=44)\n",
        "    elif link == 'probit':\n",
        "      response = tf.cast(\n",
        "          tfd.Normal(loc=linear_response, scale=scale).sample(seed=44) \u003e 0,\n",
        "                     dtype)\n",
        "    elif link == 'logit':\n",
        "      response = tfd.Bernoulli(logits=linear_response).sample(seed=44)\n",
        "    else:\n",
        "      raise ValueError('unrecognized true link: {}'.format(link))\n",
        "    return model_matrix, response, model_coefficients, mask"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "99Fk5XZKbvi4"
      },
      "source": [
        "### Note: Connect to a local runtime.\n",
        "\n",
        "In this notebook, we share data between Python and R kernels using local files.  To enable this sharing, please use runtimes on the same machine where you have permission to read and write local files."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "2EAQjTrZJqKx"
      },
      "outputs": [],
      "source": [
        "with tf.Session() as sess:\n",
        "  x_, y_, model_coefficients_true_, _ = sess.run(make_dataset(\n",
        "      n=int(1e5), d=100, link='probit'))\n",
        "\n",
        "DATA_DIR = '/tmp/glm_example'\n",
        "tf.gfile.MakeDirs(DATA_DIR)\n",
        "with tf.gfile.Open('{}/x.csv'.format(DATA_DIR), 'w') as f:\n",
        "  np.savetxt(f, x_, delimiter=',')\n",
        "with tf.gfile.Open('{}/y.csv'.format(DATA_DIR), 'w') as f:\n",
        "  np.savetxt(f, y_.astype(np.int32) + 1, delimiter=',', fmt='%d')\n",
        "with tf.gfile.Open(\n",
        "    '{}/model_coefficients_true.csv'.format(DATA_DIR), 'w') as f:\n",
        "  np.savetxt(f, model_coefficients_true_, delimiter=',')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "0P5I-aJdN6GZ"
      },
      "source": [
        "## Without L1 Regularization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VN6HfiH3bAb0"
      },
      "source": [
        "The function `tfp.glm.fit` implements Fisher scoring, which takes as some of its arguments:\n",
        "\n",
        "* `model_matrix` = $\\mathbf{x}$\n",
        "* `response` = $\\mathbf{y}$\n",
        "* `model` = callable which, given argument $\\boldsymbol{\\eta}$, returns the triple $\\left(\n",
        "{\\textbf{Mean}_T}(\\boldsymbol{\\eta}),\n",
        "{\\textbf{Var}_T}(\\boldsymbol{\\eta}),\n",
        "{\\textbf{Mean}_T}'(\\boldsymbol{\\eta})\n",
        "\\right)$.\n",
        "\n",
        "We recommend that `model` be an instance of the `tfp.glm.ExponentialFamily` class.  There are several pre-made implementations available, so for most common GLMs no custom code is necessary."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 104
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 2220,
          "status": "ok",
          "timestamp": 1538271114314,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "iXkxVBSmesjn",
        "outputId": "765d2336-cb3b-41bd-d4c9-c048de52f76f"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "is_converged: True\n",
            "    num_iter: 6\n",
            "    accuracy: 0.75356\n",
            "    deviance: -0.992467164993\n",
            "||w0-w1||_2 / (1+||w0||_2): 0.0230717559048\n"
          ]
        }
      ],
      "source": [
        "model_coefficients, linear_response, is_converged, num_iter = tfp.glm.fit(\n",
        "    model_matrix=x_,\n",
        "    response=y_,\n",
        "    model=tfp.glm.BernoulliNormalCDF())\n",
        "log_likelihood = tfp.glm.BernoulliNormalCDF().log_prob(y_, linear_response)\n",
        "\n",
        "with tf.Session() as sess:\n",
        "  [model_coefficients_, linear_response_, is_converged_, num_iter_,\n",
        "   log_likelihood_] = sess.run([\n",
        "      model_coefficients, linear_response, is_converged, num_iter,\n",
        "      log_likelihood])\n",
        "\n",
        "print(('is_converged: {}\\n'\n",
        "       '    num_iter: {}\\n'\n",
        "       '    accuracy: {}\\n'\n",
        "       '    deviance: {}\\n'\n",
        "       '||w0-w1||_2 / (1+||w0||_2): {}'\n",
        "      ).format(\n",
        "    is_converged_,\n",
        "    num_iter_,\n",
        "    np.mean((linear_response_ \u003e 0.) == y_),\n",
        "    2. * np.mean(log_likelihood_),\n",
        "    np.linalg.norm(model_coefficients_true_ - model_coefficients_, ord=2) /\n",
        "        (1. + np.linalg.norm(model_coefficients_true_, ord=2))\n",
        "    ))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "h6qexoHAJzEF"
      },
      "source": [
        "### Mathematical Details\n",
        "\n",
        "Fisher scoring is a modification of Newton's method to find the maximum-likelihood estimate\n",
        "\n",
        "$$\n",
        "\\hat\\beta\n",
        ":= \\underset{\\beta}{\\text{arg max}}\\ \\ \\ell(\\beta\\ ;\\ \\mathbf{x}, \\mathbf{y}).\n",
        "$$\n",
        "\n",
        "Vanilla Newton's method, searching for zeros of the gradient of the log-likelihood, would follow the update rule\n",
        "\n",
        "$$\n",
        "  \\beta^{(t+1)}_{\\text{Newton}}\n",
        ":=\n",
        "  \\beta^{(t)}\n",
        "  -\n",
        "  \\alpha\n",
        "  \\left(\n",
        "    \\nabla^2_\\beta\\, \\ell(\\beta\\ ;\\ \\mathbf{x}, \\mathbf{y})\n",
        "  \\right)_{\\beta = \\beta^{(t)}}^{-1}\n",
        "  \\left(\n",
        "    \\nabla_\\beta\\, \\ell(\\beta\\ ;\\ \\mathbf{x}, \\mathbf{y})\n",
        "  \\right)_{\\beta = \\beta^{(t)}}\n",
        "$$\n",
        "\n",
        "where $\\alpha \\in (0, 1]$ is a learning rate used to control the step size.\n",
        "\n",
        "In Fisher scoring, we replace the Hessian with the negative Fisher information matrix:\n",
        "\n",
        "\\begin{align*}\n",
        "  \\beta^{(t+1)}\n",
        "\u0026:=\n",
        "  \\beta^{(t)}\n",
        "  -\n",
        "  \\alpha\\,\n",
        "  \\mathbb{E}_{\n",
        "    Y_i \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x_i^\\top \\beta^{(t)}), \\phi)\n",
        "  }\n",
        "  \\left[\n",
        "    \\left(\n",
        "      \\nabla^2_\\beta\\, \\ell(\\beta\\ ;\\ \\mathbf{x}, \\mathbf{Y})\n",
        "    \\right)_{\\beta = \\beta^{(t)}}\n",
        "  \\right]^{-1}\n",
        "  \\left(\n",
        "    \\nabla_\\beta\\, \\ell(\\beta\\ ;\\ \\mathbf{x}, \\mathbf{y})\n",
        "  \\right)_{\\beta = \\beta^{(t)}} \\\\[3mm]\n",
        "\\end{align*}\n",
        "\n",
        "[Note that here $\\mathbf{Y} = (Y_i)_{i=1}^{n}$ is random, whereas $\\mathbf{y}$ is still the vector of observed responses.]\n",
        "\n",
        "By the formulas in \"Fitting GLM Parameters To Data\" below, this simplifies to\n",
        "\n",
        "\\begin{align*}\n",
        "  \\beta^{(t+1)}\n",
        "\u0026=\n",
        "  \\beta^{(t)}\n",
        "  +\n",
        "  \\alpha\n",
        "  \\left(\n",
        "    \\mathbf{x}^\\top\n",
        "    \\text{diag}\\left(\n",
        "      \\frac{\n",
        "        \\phi\\, {\\textbf{Mean}_T}'(\\mathbf{x} \\beta^{(t)})^2\n",
        "      }{\n",
        "        {\\textbf{Var}_T}(\\mathbf{x} \\beta^{(t)})\n",
        "      }\\right)\\,\n",
        "    \\mathbf{x}\n",
        "  \\right)^{-1}\n",
        "  \\left(\n",
        "    \\mathbf{x}^\\top\n",
        "    \\text{diag}\\left(\\frac{\n",
        "        {\\textbf{Mean}_T}'(\\mathbf{x} \\beta^{(t)})\n",
        "      }{\n",
        "        {\\textbf{Var}_T}(\\mathbf{x} \\beta^{(t)})\n",
        "      }\\right)\n",
        "    \\left(\\mathbf{T}(\\mathbf{y}) - {\\textbf{Mean}_T}(\\mathbf{x} \\beta^{(t)})\\right)\n",
        "  \\right).\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "076quM7tN8_1"
      },
      "source": [
        "## With L1 Regularization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "fnP3jeZOk7Y5"
      },
      "source": [
        "`tfp.glm.fit_sparse` implements a GLM fitter more suited to sparse data sets, based on the algorithm in [Yuan, Ho and Lin 2012](#1).  Its features include:\n",
        "\n",
        "* L1 regularization\n",
        "* No matrix inversions\n",
        "* Few evaluations of the gradient and Hessian.\n",
        "\n",
        "We first present an example usage of the code.  Details of the algorithm are further elaborated in \"Algorithm Details for `tfp.glm.fit_sparse`\" below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 2039
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 3299,
          "status": "ok",
          "timestamp": 1538271150420,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "v_Oky1X4ijfv",
        "outputId": "0fb5c4e2-f0f8-47f7-e423-09b90e4beef8"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "is_converged: True\n",
            "    num_iter: 1\n",
            "\n",
            "Coefficients:\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "\u003cdiv style=\"max-width:1500px;overflow:auto;\"\u003e\n",
              "\u003cstyle scoped\u003e\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "\u003c/style\u003e\n",
              "\u003ctable border=\"1\" class=\"dataframe\"\u003e\n",
              "  \u003cthead\u003e\n",
              "    \u003ctr style=\"text-align: right;\"\u003e\n",
              "      \u003cth\u003e\u003c/th\u003e\n",
              "      \u003cth\u003eLearned\u003c/th\u003e\n",
              "      \u003cth\u003eTrue\u003c/th\u003e\n",
              "    \u003c/tr\u003e\n",
              "  \u003c/thead\u003e\n",
              "  \u003ctbody\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e0\u003c/th\u003e\n",
              "      \u003ctd\u003e0.209221\u003c/td\u003e\n",
              "      \u003ctd\u003e0.220758\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e1\u003c/th\u003e\n",
              "      \u003ctd\u003e0.055613\u003c/td\u003e\n",
              "      \u003ctd\u003e0.086594\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e2\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e3\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e4\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e5\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e6\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.138448\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.153256\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e7\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e8\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.058378\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.079942\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e9\u003c/th\u003e\n",
              "      \u003ctd\u003e0.085470\u003c/td\u003e\n",
              "      \u003ctd\u003e0.109012\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e10\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.113677\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.137061\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e11\u003c/th\u003e\n",
              "      \u003ctd\u003e0.163041\u003c/td\u003e\n",
              "      \u003ctd\u003e0.174420\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e12\u003c/th\u003e\n",
              "      \u003ctd\u003e0.141059\u003c/td\u003e\n",
              "      \u003ctd\u003e0.158109\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e13\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e14\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.246402\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.242609\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e15\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.107397\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.123367\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e16\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.002965\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e17\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.037789\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.067560\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e18\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e19\u003c/th\u003e\n",
              "      \u003ctd\u003e0.228254\u003c/td\u003e\n",
              "      \u003ctd\u003e0.230469\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e20\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e21\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e22\u003c/th\u003e\n",
              "      \u003ctd\u003e0.024156\u003c/td\u003e\n",
              "      \u003ctd\u003e0.049624\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e23\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e24\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e25\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.109751\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.128736\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e26\u003c/th\u003e\n",
              "      \u003ctd\u003e0.062569\u003c/td\u003e\n",
              "      \u003ctd\u003e0.085314\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e27\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.235091\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.240482\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e28\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e29\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e...\u003c/th\u003e\n",
              "      \u003ctd\u003e...\u003c/td\u003e\n",
              "      \u003ctd\u003e...\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e70\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e71\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e72\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.100460\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.117774\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e73\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.134124\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.149282\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e74\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e75\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e76\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.181901\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.195186\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e77\u003c/th\u003e\n",
              "      \u003ctd\u003e0.222551\u003c/td\u003e\n",
              "      \u003ctd\u003e0.229656\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e78\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e79\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.229750\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.225494\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e80\u003c/th\u003e\n",
              "      \u003ctd\u003e0.196564\u003c/td\u003e\n",
              "      \u003ctd\u003e0.200442\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e81\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e82\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.029118\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.053469\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e83\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e84\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e85\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e86\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e87\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e88\u003c/th\u003e\n",
              "      \u003ctd\u003e0.155175\u003c/td\u003e\n",
              "      \u003ctd\u003e0.171693\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e89\u003c/th\u003e\n",
              "      \u003ctd\u003e0.200711\u003c/td\u003e\n",
              "      \u003ctd\u003e0.206408\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e90\u003c/th\u003e\n",
              "      \u003ctd\u003e0.039990\u003c/td\u003e\n",
              "      \u003ctd\u003e0.063434\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e91\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e92\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e93\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e94\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e95\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e96\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.014878\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.046784\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e97\u003c/th\u003e\n",
              "      \u003ctd\u003e0.055308\u003c/td\u003e\n",
              "      \u003ctd\u003e0.081945\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e98\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.076264\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.104929\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e99\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.138370\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.153871\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "  \u003c/tbody\u003e\n",
              "\u003c/table\u003e\n",
              "\u003cp\u003e100 rows × 2 columns\u003c/p\u003e\n",
              "\u003c/div\u003e"
            ],
            "text/plain": [
              "     Learned      True\n",
              "0   0.209221  0.220758\n",
              "1   0.055613  0.086594\n",
              "2   0.000000  0.000000\n",
              "3   0.000000  0.000000\n",
              "4   0.000000  0.000000\n",
              "5   0.000000  0.000000\n",
              "6  -0.138448 -0.153256\n",
              "7   0.000000  0.000000\n",
              "8  -0.058378 -0.079942\n",
              "9   0.085470  0.109012\n",
              "10 -0.113677 -0.137061\n",
              "11  0.163041  0.174420\n",
              "12  0.141059  0.158109\n",
              "13  0.000000  0.000000\n",
              "14 -0.246402 -0.242609\n",
              "15 -0.107397 -0.123367\n",
              "16  0.000000  0.002965\n",
              "17 -0.037789 -0.067560\n",
              "18  0.000000  0.000000\n",
              "19  0.228254  0.230469\n",
              "20  0.000000  0.000000\n",
              "21  0.000000  0.000000\n",
              "22  0.024156  0.049624\n",
              "23  0.000000  0.000000\n",
              "24  0.000000  0.000000\n",
              "25 -0.109751 -0.128736\n",
              "26  0.062569  0.085314\n",
              "27 -0.235091 -0.240482\n",
              "28  0.000000  0.000000\n",
              "29  0.000000  0.000000\n",
              "..       ...       ...\n",
              "70  0.000000  0.000000\n",
              "71  0.000000  0.000000\n",
              "72 -0.100460 -0.117774\n",
              "73 -0.134124 -0.149282\n",
              "74  0.000000  0.000000\n",
              "75  0.000000  0.000000\n",
              "76 -0.181901 -0.195186\n",
              "77  0.222551  0.229656\n",
              "78  0.000000  0.000000\n",
              "79 -0.229750 -0.225494\n",
              "80  0.196564  0.200442\n",
              "81  0.000000  0.000000\n",
              "82 -0.029118 -0.053469\n",
              "83  0.000000  0.000000\n",
              "84  0.000000  0.000000\n",
              "85  0.000000  0.000000\n",
              "86  0.000000  0.000000\n",
              "87  0.000000  0.000000\n",
              "88  0.155175  0.171693\n",
              "89  0.200711  0.206408\n",
              "90  0.039990  0.063434\n",
              "91  0.000000  0.000000\n",
              "92  0.000000  0.000000\n",
              "93  0.000000  0.000000\n",
              "94  0.000000  0.000000\n",
              "95  0.000000  0.000000\n",
              "96 -0.014878 -0.046784\n",
              "97  0.055308  0.081945\n",
              "98 -0.076264 -0.104929\n",
              "99 -0.138370 -0.153871\n",
              "\n",
              "[100 rows x 2 columns]"
            ]
          },
          "execution_count": 7,
          "metadata": {
            "tags": []
          },
          "output_type": "execute_result"
        }
      ],
      "source": [
        "with tf.Session() as sess:\n",
        "  model = tfp.glm.Bernoulli()\n",
        "  model_coefficients_start = tf.zeros(x_.shape[-1], np.float32)\n",
        "\n",
        "  with tf.variable_scope(tf.get_variable_scope(), reuse=tf.AUTO_REUSE):\n",
        "    model_coefficients, is_converged, num_iter = tfp.glm.fit_sparse(\n",
        "        model_matrix=tf.convert_to_tensor(x_),\n",
        "        response=tf.convert_to_tensor(y_),\n",
        "        model=model,\n",
        "        model_coefficients_start=model_coefficients_start,\n",
        "        l1_regularizer=800.,\n",
        "        l2_regularizer=None,\n",
        "        maximum_iterations=10,\n",
        "        maximum_full_sweeps_per_iteration=10,\n",
        "        tolerance=1e-6,\n",
        "        learning_rate=None)\n",
        "\n",
        "  init_op = tf.global_variables_initializer()\n",
        "  sess.run([init_op])\n",
        "  model_coefficients_, is_converged_, num_iter_ = sess.run([\n",
        "      model_coefficients, is_converged, num_iter])\n",
        "  coefs_comparison = pd.DataFrame({\n",
        "    'Learned': model_coefficients_,\n",
        "    'True': model_coefficients_true_,\n",
        "  })\n",
        "  \n",
        "print(('is_converged: {}\\n'\n",
        "       '    num_iter: {}\\n\\n'\n",
        "       'Coefficients:').format(\n",
        "    is_converged_,\n",
        "    num_iter_))\n",
        "coefs_comparison"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "DrJC2J1YbR5L"
      },
      "source": [
        "Note that the learned coefficients have the same sparsity pattern as the true coefficients."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "hQ7SzrPZMpke"
      },
      "outputs": [],
      "source": [
        "# Save the learned coefficients to a file.\n",
        "with tf.gfile.Open('{}/model_coefficients_prox.csv'.format(DATA_DIR), 'w') as f:\n",
        "  np.savetxt(f, model_coefficients_, delimiter=',')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VW9NgB1Zisqh"
      },
      "source": [
        "### Compare to R's `glmnet`\n",
        "\n",
        "We compare the output of coordinatewise proximal gradient descent to that of R's `glmnet`, which uses a similar algorithm."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Aptz7SWwkd5v"
      },
      "source": [
        "#### NOTE: To execute this section, you must switch to an R colab runtime."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "RS1H3n53h9qc"
      },
      "outputs": [],
      "source": [
        "suppressMessages({\n",
        "  library('glmnet')\n",
        "})"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "2X6zKSaxie7I"
      },
      "outputs": [],
      "source": [
        "data_dir \u003c- '/tmp/glm_example'\n",
        "x \u003c- as.matrix(read.csv(paste(data_dir, '/x.csv', sep=''),\n",
        "                        header=FALSE))\n",
        "y \u003c- as.matrix(read.csv(paste(data_dir, '/y.csv', sep=''),\n",
        "                        header=FALSE, colClasses='integer'))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "Eb31LbhRjsSz"
      },
      "outputs": [],
      "source": [
        "fit \u003c- glmnet(\n",
        "x = x,\n",
        "y = y,\n",
        "family = \"binomial\",  # Logistic regression\n",
        "alpha = 1,  # corresponds to l1_weight = 1, l2_weight = 0\n",
        "standardize = FALSE,\n",
        "intercept = FALSE,\n",
        "thresh = 1e-30,\n",
        "type.logistic = \"Newton\"\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "HTN4RKQbhlCm"
      },
      "outputs": [],
      "source": [
        "write.csv(as.matrix(coef(fit, 0.008)),\n",
        "          paste(data_dir, '/model_coefficients_glmnet.csv', sep=''),\n",
        "          row.names=FALSE)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "vsrEKgUGjGjf"
      },
      "source": [
        "#### Compare R, TFP and true coefficients (Note: back to Python kernel)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {},
        "colab_type": "code",
        "id": "lCOlGo_4i2sb"
      },
      "outputs": [],
      "source": [
        "DATA_DIR = '/tmp/glm_example'\n",
        "with tf.gfile.Open('{}/model_coefficients_glmnet.csv'.format(DATA_DIR),\n",
        "                   'r') as f:\n",
        "  model_coefficients_glmnet_ = np.loadtxt(f,\n",
        "                                   skiprows=2  # Skip column name and intercept\n",
        "                               )\n",
        "\n",
        "with tf.gfile.Open('{}/model_coefficients_prox.csv'.format(DATA_DIR), 'r') as f:\n",
        "  model_coefficients_prox_ = np.loadtxt(f)\n",
        "\n",
        "with tf.gfile.Open(\n",
        "    '{}/model_coefficients_true.csv'.format(DATA_DIR), 'r') as f:\n",
        "  model_coefficients_true_ = np.loadtxt(f)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 1969
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 218,
          "status": "ok",
          "timestamp": 1538271869238,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "4l-SZ85lnKg5",
        "outputId": "e53b8078-2338-4168-c54d-18975e9a26a7"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "\u003cdiv style=\"max-width:1500px;overflow:auto;\"\u003e\n",
              "\u003cstyle scoped\u003e\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "\u003c/style\u003e\n",
              "\u003ctable border=\"1\" class=\"dataframe\"\u003e\n",
              "  \u003cthead\u003e\n",
              "    \u003ctr style=\"text-align: right;\"\u003e\n",
              "      \u003cth\u003e\u003c/th\u003e\n",
              "      \u003cth\u003eR\u003c/th\u003e\n",
              "      \u003cth\u003eTFP\u003c/th\u003e\n",
              "      \u003cth\u003eTrue\u003c/th\u003e\n",
              "    \u003c/tr\u003e\n",
              "  \u003c/thead\u003e\n",
              "  \u003ctbody\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e0\u003c/th\u003e\n",
              "      \u003ctd\u003e0.271080\u003c/td\u003e\n",
              "      \u003ctd\u003e0.209221\u003c/td\u003e\n",
              "      \u003ctd\u003e0.220758\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e1\u003c/th\u003e\n",
              "      \u003ctd\u003e0.073014\u003c/td\u003e\n",
              "      \u003ctd\u003e0.055613\u003c/td\u003e\n",
              "      \u003ctd\u003e0.086594\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e2\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e3\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e4\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e5\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e6\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.178551\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.138448\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.153256\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e7\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e8\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.075707\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.058378\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.079942\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e9\u003c/th\u003e\n",
              "      \u003ctd\u003e0.110585\u003c/td\u003e\n",
              "      \u003ctd\u003e0.085470\u003c/td\u003e\n",
              "      \u003ctd\u003e0.109012\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e10\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.146996\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.113677\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.137061\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e11\u003c/th\u003e\n",
              "      \u003ctd\u003e0.211116\u003c/td\u003e\n",
              "      \u003ctd\u003e0.163041\u003c/td\u003e\n",
              "      \u003ctd\u003e0.174420\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e12\u003c/th\u003e\n",
              "      \u003ctd\u003e0.182786\u003c/td\u003e\n",
              "      \u003ctd\u003e0.141059\u003c/td\u003e\n",
              "      \u003ctd\u003e0.158109\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e13\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e14\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.319566\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.246402\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.242609\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e15\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.139507\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.107397\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.123367\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e16\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.002965\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e17\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.048688\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.037789\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.067560\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e18\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e19\u003c/th\u003e\n",
              "      \u003ctd\u003e0.296542\u003c/td\u003e\n",
              "      \u003ctd\u003e0.228254\u003c/td\u003e\n",
              "      \u003ctd\u003e0.230469\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e20\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e21\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e22\u003c/th\u003e\n",
              "      \u003ctd\u003e0.031309\u003c/td\u003e\n",
              "      \u003ctd\u003e0.024156\u003c/td\u003e\n",
              "      \u003ctd\u003e0.049624\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e23\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e24\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e25\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.141382\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.109751\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.128736\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e26\u003c/th\u003e\n",
              "      \u003ctd\u003e0.081012\u003c/td\u003e\n",
              "      \u003ctd\u003e0.062569\u003c/td\u003e\n",
              "      \u003ctd\u003e0.085314\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e27\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.304464\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.235091\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.240482\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e28\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e29\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e...\u003c/th\u003e\n",
              "      \u003ctd\u003e...\u003c/td\u003e\n",
              "      \u003ctd\u003e...\u003c/td\u003e\n",
              "      \u003ctd\u003e...\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e70\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e71\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e72\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.130267\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.100460\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.117774\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e73\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.173815\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.134124\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.149282\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e74\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e75\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e76\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.235598\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.181901\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.195186\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e77\u003c/th\u003e\n",
              "      \u003ctd\u003e0.288192\u003c/td\u003e\n",
              "      \u003ctd\u003e0.222551\u003c/td\u003e\n",
              "      \u003ctd\u003e0.229656\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e78\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e79\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.295628\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.229750\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.225494\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e80\u003c/th\u003e\n",
              "      \u003ctd\u003e0.254703\u003c/td\u003e\n",
              "      \u003ctd\u003e0.196564\u003c/td\u003e\n",
              "      \u003ctd\u003e0.200442\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e81\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e82\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.037550\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.029118\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.053469\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e83\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e84\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e85\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e86\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e87\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e88\u003c/th\u003e\n",
              "      \u003ctd\u003e0.201640\u003c/td\u003e\n",
              "      \u003ctd\u003e0.155175\u003c/td\u003e\n",
              "      \u003ctd\u003e0.171693\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e89\u003c/th\u003e\n",
              "      \u003ctd\u003e0.259279\u003c/td\u003e\n",
              "      \u003ctd\u003e0.200711\u003c/td\u003e\n",
              "      \u003ctd\u003e0.206408\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e90\u003c/th\u003e\n",
              "      \u003ctd\u003e0.052837\u003c/td\u003e\n",
              "      \u003ctd\u003e0.039990\u003c/td\u003e\n",
              "      \u003ctd\u003e0.063434\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e91\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e92\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e93\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e94\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e95\u003c/th\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "      \u003ctd\u003e0.000000\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e96\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.019880\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.014878\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.046784\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e97\u003c/th\u003e\n",
              "      \u003ctd\u003e0.071567\u003c/td\u003e\n",
              "      \u003ctd\u003e0.055308\u003c/td\u003e\n",
              "      \u003ctd\u003e0.081945\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e98\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.098354\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.076264\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.104929\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "    \u003ctr\u003e\n",
              "      \u003cth\u003e99\u003c/th\u003e\n",
              "      \u003ctd\u003e-0.179619\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.138370\u003c/td\u003e\n",
              "      \u003ctd\u003e-0.153871\u003c/td\u003e\n",
              "    \u003c/tr\u003e\n",
              "  \u003c/tbody\u003e\n",
              "\u003c/table\u003e\n",
              "\u003cp\u003e100 rows × 3 columns\u003c/p\u003e\n",
              "\u003c/div\u003e"
            ],
            "text/plain": [
              "           R       TFP      True\n",
              "0   0.271080  0.209221  0.220758\n",
              "1   0.073014  0.055613  0.086594\n",
              "2   0.000000  0.000000  0.000000\n",
              "3   0.000000  0.000000  0.000000\n",
              "4   0.000000  0.000000  0.000000\n",
              "5   0.000000  0.000000  0.000000\n",
              "6  -0.178551 -0.138448 -0.153256\n",
              "7   0.000000  0.000000  0.000000\n",
              "8  -0.075707 -0.058378 -0.079942\n",
              "9   0.110585  0.085470  0.109012\n",
              "10 -0.146996 -0.113677 -0.137061\n",
              "11  0.211116  0.163041  0.174420\n",
              "12  0.182786  0.141059  0.158109\n",
              "13  0.000000  0.000000  0.000000\n",
              "14 -0.319566 -0.246402 -0.242609\n",
              "15 -0.139507 -0.107397 -0.123367\n",
              "16  0.000000  0.000000  0.002965\n",
              "17 -0.048688 -0.037789 -0.067560\n",
              "18  0.000000  0.000000  0.000000\n",
              "19  0.296542  0.228254  0.230469\n",
              "20  0.000000  0.000000  0.000000\n",
              "21  0.000000  0.000000  0.000000\n",
              "22  0.031309  0.024156  0.049624\n",
              "23  0.000000  0.000000  0.000000\n",
              "24  0.000000  0.000000  0.000000\n",
              "25 -0.141382 -0.109751 -0.128736\n",
              "26  0.081012  0.062569  0.085314\n",
              "27 -0.304464 -0.235091 -0.240482\n",
              "28  0.000000  0.000000  0.000000\n",
              "29  0.000000  0.000000  0.000000\n",
              "..       ...       ...       ...\n",
              "70  0.000000  0.000000  0.000000\n",
              "71  0.000000  0.000000  0.000000\n",
              "72 -0.130267 -0.100460 -0.117774\n",
              "73 -0.173815 -0.134124 -0.149282\n",
              "74  0.000000  0.000000  0.000000\n",
              "75  0.000000  0.000000  0.000000\n",
              "76 -0.235598 -0.181901 -0.195186\n",
              "77  0.288192  0.222551  0.229656\n",
              "78  0.000000  0.000000  0.000000\n",
              "79 -0.295628 -0.229750 -0.225494\n",
              "80  0.254703  0.196564  0.200442\n",
              "81  0.000000  0.000000  0.000000\n",
              "82 -0.037550 -0.029118 -0.053469\n",
              "83  0.000000  0.000000  0.000000\n",
              "84  0.000000  0.000000  0.000000\n",
              "85  0.000000  0.000000  0.000000\n",
              "86  0.000000  0.000000  0.000000\n",
              "87  0.000000  0.000000  0.000000\n",
              "88  0.201640  0.155175  0.171693\n",
              "89  0.259279  0.200711  0.206408\n",
              "90  0.052837  0.039990  0.063434\n",
              "91  0.000000  0.000000  0.000000\n",
              "92  0.000000  0.000000  0.000000\n",
              "93  0.000000  0.000000  0.000000\n",
              "94  0.000000  0.000000  0.000000\n",
              "95  0.000000  0.000000  0.000000\n",
              "96 -0.019880 -0.014878 -0.046784\n",
              "97  0.071567  0.055308  0.081945\n",
              "98 -0.098354 -0.076264 -0.104929\n",
              "99 -0.179619 -0.138370 -0.153871\n",
              "\n",
              "[100 rows x 3 columns]"
            ]
          },
          "execution_count": 10,
          "metadata": {
            "tags": []
          },
          "output_type": "execute_result"
        }
      ],
      "source": [
        "coefs_comparison = pd.DataFrame({\n",
        "    'TFP': model_coefficients_prox_,\n",
        "    'R': model_coefficients_glmnet_,\n",
        "    'True': model_coefficients_true_,\n",
        "})\n",
        "coefs_comparison"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Rfv0GVXqY74Y"
      },
      "source": [
        "# Algorithm Details for `tfp.glm.fit_sparse`\n",
        "\n",
        "We present the algorithm as a sequence of three modifications to Newton's method.  In each one, the update rule for $\\beta$ is based on a vector $s$ and a matrix $H$ which approximate the gradient and Hessian of the log-likelihood.  In step $t$, we choose a coordinate $j^{(t)}$ to change, and we update $\\beta$ according to the update rule:\n",
        "\n",
        "\\begin{align*}\n",
        "  u^{(t)}\n",
        "\u0026:=\n",
        "  \\frac{\n",
        "    \\left(\n",
        "      s^{(t)}\n",
        "    \\right)_{j^{(t)}}\n",
        "  }{\n",
        "    \\left(\n",
        "      H^{(t)}\n",
        "    \\right)_{j^{(t)},\\, j^{(t)}}\n",
        "  }\n",
        "\\\\[3mm]\n",
        "  \\beta^{(t+1)}\n",
        "\u0026:=\n",
        "  \\beta^{(t)}\n",
        "  -\n",
        "  \\alpha\\,\n",
        "  u^{(t)}\n",
        "  \\,\\text{onehot}(j^{(t)})\n",
        "\\end{align*}\n",
        "\n",
        "This update is a Newton-like step with learning rate $\\alpha$.  Except for the final piece (L1 regularization), the modifications below differ only in how they update $s$ and $H$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "fH7C1xBWUV7_"
      },
      "source": [
        "## Starting point: Coordinatewise Newton's method\n",
        "\n",
        "In coordinatewise Newton's method, we set $s$ and $H$ to the true gradient and Hessian of the log-likelihood:\n",
        "\n",
        "\\begin{align*}\n",
        "  s^{(t)}_{\\text{vanilla}}\n",
        "\u0026:=\n",
        "  \\left(\n",
        "    \\nabla_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "  \\right)_{\\beta = \\beta^{(t)}}\n",
        "\\\\\n",
        "  H^{(t)}_{\\text{vanilla}}\n",
        "\u0026:=\n",
        "  \\left(\n",
        "    \\nabla^2_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "  \\right)_{\\beta = \\beta^{(t)}}\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "6rJZD6iyUl0v"
      },
      "source": [
        "## Fewer evaluations of the gradient and Hessian\n",
        "\n",
        "The gradient and Hessian of the log-likelihood are often expensive to compute, so it is often worthwhile to approximate them.  We can do so as follows:\n",
        "\n",
        "* Usually, approximate the Hessian as locally constant and approximate the gradient to first order using the (approximate) Hessian:\n",
        "\n",
        "\\begin{align*}\n",
        "  H_{\\text{approx}}^{(t+1)}\n",
        "\u0026:=\n",
        "  H^{(t)}\n",
        "\\\\\n",
        "  s_{\\text{approx}}^{(t+1)}\n",
        "\u0026:=\n",
        "  s^{(t)}\n",
        "  +\n",
        "  H^{(t)}\n",
        "  \\left(\n",
        "    \\beta^{(t+1)} - \\beta^{(t)}\n",
        "  \\right)\n",
        "\\end{align*}\n",
        "\n",
        "* Occasionally, perform a \"vanilla\" update step as above, setting $s^{(t+1)}$ to the exact gradient and $H^{(t+1)}$ to the exact Hessian of the log-likelihood, evaluated at $\\beta^{(t+1)}$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "rfvvyaVnUqIQ"
      },
      "source": [
        "## Substitute negative Fisher information for Hessian\n",
        "\n",
        "To further reduce the cost of the vanilla update steps, we can set $H$ to the negative Fisher information matrix (efficiently computable using the formulas in \"Fitting GLM Parameters to Data\" below) rather than the exact Hessian:\n",
        "\n",
        "\\begin{align*}\n",
        "  H_{\\text{Fisher}}^{(t+1)}\n",
        "\u0026:=\n",
        "  \\mathbb{E}_{Y_i \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x_i^\\top \\beta^{(t+1)}), \\phi)}\n",
        "  \\left[\n",
        "    \\left(\n",
        "      \\nabla_\\beta^2\\, \\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{Y})\n",
        "    \\right)_{\\beta = \\beta^{(t+1)}}\n",
        "  \\right]\n",
        "\\\\\n",
        "\u0026=\n",
        "  -\\mathbf{x}^\\top\n",
        "  \\,\\text{diag}\\left(\n",
        "    \\frac{\n",
        "      \\phi\\, {\\textbf{Mean}_T}'(\\mathbf{x} \\beta^{(t+1)})^2\n",
        "    }{\n",
        "      {\\textbf{Var}_T}(\\mathbf{x} \\beta^{(t+1)})\n",
        "    }\\right)\\,\n",
        "  \\mathbf{x}\n",
        "\\\\\n",
        "  s_{\\text{Fisher}}^{(t+1)}\n",
        "\u0026:=\n",
        "  s_{\\text{vanilla}}^{(t+1)}\n",
        "\\\\\n",
        "\u0026=\n",
        "  \\left(\n",
        "    \\mathbf{x}^\\top\n",
        "    \\,\\text{diag}\\left(\\frac{\n",
        "        {\\textbf{Mean}_T}'(\\mathbf{x} \\beta^{(t+1)})\n",
        "      }{\n",
        "        {\\textbf{Var}_T}(\\mathbf{x} \\beta^{(t+1)})\n",
        "      }\\right)\n",
        "    \\left(\\mathbf{T}(\\mathbf{y}) - {\\textbf{Mean}_T}(\\mathbf{x} \\beta^{(t+1)})\\right)\n",
        "  \\right)\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "DTH07xYpWGcR"
      },
      "source": [
        "## L1 Regularization via Proximal Gradient Descent\n",
        "\n",
        "To incorporate L1 regularization, we replace the update rule\n",
        "\n",
        "$$\n",
        "  \\beta^{(t+1)}\n",
        ":=\n",
        "  \\beta^{(t)}\n",
        "  -\n",
        "  \\alpha\\,\n",
        "  u^{(t)}\n",
        "  \\,\\text{onehot}(j^{(t)})\n",
        "$$\n",
        "\n",
        "with the more general update rule\n",
        "\n",
        "\\begin{align*}\n",
        "  \\gamma^{(t)}\n",
        "\u0026:=\n",
        "  -\\frac{\\alpha\\, r_{\\text{L1}}}{\\left(H^{(t)}\\right)_{j^{(t)},\\, j^{(t)}}}\n",
        "\\\\[2mm]\n",
        "  \\left(\\beta_{\\text{reg}}^{(t+1)}\\right)_j\n",
        "\u0026:=\n",
        "  \\begin{cases}\n",
        "    \\beta^{(t+1)}_j\n",
        "      \u0026\\text{if } j \\neq j^{(t)} \\\\\n",
        "    \\text{SoftThreshold} \\left(\n",
        "      \\beta^{(t)}_j - \\alpha\\, u^{(t)}\n",
        "      ,\\ \n",
        "      \\gamma^{(t)}\n",
        "    \\right)\n",
        "      \u0026\\text{if } j = j^{(t)}\n",
        "  \\end{cases}\n",
        "\\end{align*}\n",
        "\n",
        "where $r_{\\text{L1}} \u003e 0$ is a supplied constant (the L1 regularization coefficient) and $\\text{SoftThreshold}$ is the soft thresholding operator, defined by\n",
        "\n",
        "$$\n",
        "\\text{SoftThreshold}(\\beta, \\gamma)\n",
        ":=\n",
        "\\begin{cases}\n",
        "\\beta + \\gamma\n",
        "  \u0026\\text{if } \\beta \u003c -\\gamma\n",
        "\\\\\n",
        "0\n",
        "  \u0026\\text{if } -\\gamma \\leq \\beta \\leq \\gamma\n",
        "\\\\\n",
        "\\beta - \\gamma\n",
        "  \u0026\\text{if } \\beta \u003e \\gamma.\n",
        "\\end{cases}\n",
        "$$\n",
        "\n",
        "This update rule has the following two inspirational properties, which we explain below:\n",
        "\n",
        "1. In the limiting case $r_{\\text{L1}} \\to 0$ (i.e., no L1 regularization), this update rule is identical to the original update rule.\n",
        "\n",
        "1. This update rule can be interpreted as applying a proximity operator whose fixed point is the solution to the L1-regularized minimization problem\n",
        "\n",
        "$$\n",
        "\\underset{\\beta - \\beta^{(t)} \\in \\text{span}\\{ \\text{onehot}(j^{(t)}) \\}}{\\text{arg min}}\n",
        "\\left(\n",
        "  -\\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "  + r_{\\text{L1}} \\left\\lVert \\beta \\right\\rVert_1\n",
        "\\right).\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "CSs7_osNPLVt"
      },
      "source": [
        "### Degenerate case $r_{\\text{L1}} = 0$ recovers the original update rule\n",
        "\n",
        "To see (1), note that if $r_{\\text{L1}} = 0$ then $\\gamma^{(t)} = 0$, hence\n",
        "\n",
        "\\begin{align*}\n",
        "  \\left(\\beta_{\\text{reg}}^{(t+1)}\\right)_{j^{(t)}}\n",
        "\u0026=\n",
        "  \\text{SoftThreshold} \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}} - \\alpha\\, u^{(t)}\n",
        "    ,\\ \n",
        "    0\n",
        "  \\right)\n",
        "\\\\\n",
        "\u0026=\n",
        "  \\beta^{(t)}_{j^{(t)}} - \\alpha\\, u^{(t)}.\n",
        "\\end{align*}\n",
        "\n",
        "Hence\n",
        "\n",
        "\\begin{align*}\n",
        "  \\beta_{\\text{reg}}^{(t+1)}\n",
        "\u0026=\n",
        "  \\beta^{(t)} - \\alpha\\, u^{(t)} \\,\\text{onehot}(j^{(t)})\n",
        "\\\\\n",
        "\u0026=\n",
        "  \\beta^{(t+1)}.\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "EiHy_0NIPT5f"
      },
      "source": [
        "### Proximity operator whose fixed point is the regularized MLE\n",
        "\n",
        "To see (2), first note (see [Wikipedia](#3)) that for any $\\gamma \u003e 0$, the update rule\n",
        "\n",
        "$$\n",
        "  \\left(\\beta_{\\text{exact-prox}, \\gamma}^{(t+1)}\\right)_{j^{(t)}}\n",
        ":=\n",
        "  \\text{prox}_{\\gamma \\lVert \\cdot \\rVert_1}\n",
        "  \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}} + \\frac{\\gamma}{r_{\\text{L1}}}\n",
        "    \\left(\n",
        "      \\left(\n",
        "        \\nabla_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "      \\right)_{\\beta = \\beta^{(t)}}\n",
        "    \\right)_{j^{(t)}}\n",
        "  \\right)\n",
        "$$\n",
        "\n",
        "satisfies (2), where $\\text{prox}$ is the proximity operator (see [Yu](#4), where this operator is denoted $\\mathsf{P}$).  The right-hand side of the above equation is computed [here](#2):\n",
        "\n",
        "$$\n",
        "  \\left(\\beta_{\\text{exact-prox}, \\gamma}^{(t+1)}\\right)_{j^{(t)}}\n",
        "=\n",
        "  \\text{SoftThreshold} \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}}\n",
        "    +\n",
        "    \\frac{\\gamma}{r_{\\text{L1}}}\n",
        "    \\left(\n",
        "      \\left(\n",
        "        \\nabla_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "      \\right)_{\\beta = \\beta^{(t)}}\n",
        "    \\right)_{j^{(t)}}\n",
        "    ,\\ \n",
        "    \\gamma\n",
        "  \\right).\n",
        "$$\n",
        "\n",
        "In particular, setting\n",
        "$\\gamma = \\gamma^{(t)} = -\\frac{\\alpha\\, r_{\\text{L1}}}{\\left(H^{(t)}\\right)_{j^{(t)}, j^{(t)}}}$\n",
        "(note that $\\gamma^{(t)} \u003e 0$ as long as the negative log-likelihood is convex),\n",
        "we obtain the update rule\n",
        "\n",
        "$$\n",
        "  \\left(\\beta_{\\text{exact-prox}, \\gamma^{(t)}}^{(t+1)}\\right)_{j^{(t)}}\n",
        "=\n",
        "  \\text{SoftThreshold} \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}}\n",
        "    -\n",
        "    \\alpha\n",
        "    \\frac{\n",
        "      \\left(\n",
        "        \\left(\n",
        "          \\nabla_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "        \\right)_{\\beta = \\beta^{(t)}}\n",
        "      \\right)_{j^{(t)}}\n",
        "    }{\n",
        "      \\left(H^{(t)}\\right)_{j^{(t)}, j^{(t)}}\n",
        "    }\n",
        "    ,\\ \n",
        "    \\gamma^{(t)}\n",
        "  \\right).\n",
        "$$\n",
        "\n",
        "We then replace the exact gradient\n",
        "$\\left(\n",
        "        \\nabla_\\beta\\, \\ell(\\beta \\,;\\, \\mathbf{x}, \\mathbf{y})\n",
        "      \\right)_{\\beta = \\beta^{(t)}}$\n",
        "with its approximation $s^{(t)}$, obtaining\n",
        "\n",
        "\\begin{align*}\n",
        "  \\left(\\beta_{\\text{exact-prox}, \\gamma^{(t)}}^{(t+1)}\\right)_{j^{(t)}}\n",
        "\u0026\\approx\n",
        "  \\text{SoftThreshold} \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}}\n",
        "    -\n",
        "    \\alpha\n",
        "    \\frac{\n",
        "      \\left(s^{(t)}\\right)_{j^{(t)}}\n",
        "    }{\n",
        "      \\left(H^{(t)}\\right)_{j^{(t)}, j^{(t)}}\n",
        "    }\n",
        "    ,\\ \n",
        "    \\gamma^{(t)}\n",
        "  \\right)\n",
        "\\\\\n",
        "\u0026=\n",
        "    \\text{SoftThreshold} \\left(\n",
        "    \\beta^{(t)}_{j^{(t)}}\n",
        "    -\n",
        "    \\alpha\\,\n",
        "    u^{(t)}\n",
        "    ,\\ \n",
        "    \\gamma^{(t)}\n",
        "  \\right).\n",
        "\\end{align*}\n",
        "\n",
        "Hence\n",
        "\n",
        "$$\n",
        "  \\beta_{\\text{exact-prox}, \\gamma^{(t)}}^{(t+1)}\n",
        "\\approx\n",
        "  \\beta_{\\text{reg}}^{(t+1)}.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "P7YOOrmI8j0L"
      },
      "source": [
        "# Derivation of GLM Facts\n",
        "\n",
        "In this section we state in full detail and derive the results about GLMs that are used in the preceding sections.  Then, we use TensorFlow's `gradients` to numerically verify the derived formulas for gradient of the log-likelihood and Fisher information."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "lkHZyhuAIW-p"
      },
      "source": [
        "## Score and Fisher information"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "bbyYy0bE8pOK"
      },
      "source": [
        "Consider a family of probability distributions parameterized by parameter vector $\\theta$, having probability densities $\\left\\{p(\\cdot | \\theta)\\right\\}_{\\theta \\in \\mathcal{T}}$.  The **score** of an outcome $y$ at parameter vector $\\theta_0$ is defined to be the gradient of the log likelihood of $y$ (evaluated at $\\theta_0$), that is,\n",
        "\n",
        "$$\n",
        "\\text{score}(y, \\theta_0) := \\left[\\nabla_\\theta\\, \\log p(y | \\theta)\\right]_{\\theta=\\theta_0}.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "IYGaMPIx8uOc"
      },
      "source": [
        "### Claim: Expectation of the score is zero\n",
        "Under mild regularity conditions (permitting us to pass differentiation under the integral),\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\\text{score}(Y, \\theta_0)\\right] = 0.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "b3H-wNmJ800R"
      },
      "source": [
        "#### Proof\n",
        "We have\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\\text{score}(Y, \\theta_0)\\right]\n",
        "\u0026:=\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\\left(\\nabla_\\theta \\log p(Y|\\theta)\\right)_{\\theta=\\theta_0}\\right] \\\\\n",
        "\u0026\\stackrel{\\text{(1)}}{=} \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\\frac{\\left(\\nabla_\\theta p(Y|\\theta)\\right)_{\\theta=\\theta_0}}{p(Y|\\theta=\\theta_0)}\\right] \\\\\n",
        "\u0026\\stackrel{\\text{(2)}}{=} \\int_{\\mathcal{Y}} \\left[\\frac{\\left(\\nabla_\\theta p(y|\\theta)\\right)_{\\theta=\\theta_0}}{p(y|\\theta=\\theta_0)}\\right] p(y | \\theta=\\theta_0)\\, dy \\\\\n",
        "\u0026= \\int_{\\mathcal{Y}} \\left(\\nabla_\\theta p(y|\\theta)\\right)_{\\theta=\\theta_0}\\, dy \\\\\n",
        "\u0026\\stackrel{\\text{(3)}}{=} \\left[\\nabla_\\theta \\left(\\int_{\\mathcal{Y}} p(y|\\theta)\\, dy\\right) \\right]_{\\theta=\\theta_0} \\\\\n",
        "\u0026\\stackrel{\\text{(4)}}{=} \\left[\\nabla_\\theta\\, 1 \\right]_{\\theta=\\theta_0} \\\\\n",
        "\u0026= 0,\n",
        "\\end{align*}\n",
        "\n",
        "where we have used: (1) chain rule for differentiation, (2) definition of expectation, (3) passing differentiation under the integral sign (using the regularity conditions), (4) the integral of a probability density is 1."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "1Y1DPVOI9OT2"
      },
      "source": [
        "### Claim (Fisher information): Variance of the score equals negative expected Hessian of the log likelihood\n",
        "\n",
        "Under mild regularity conditions (permitting us to pass differentiation under the integral),\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "\\text{score}(Y, \\theta_0) \\text{score}(Y, \\theta_0)^\\top\n",
        "\\right]\n",
        "=\n",
        "-\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "\\left(\\nabla_\\theta^2 \\log p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "\\right]\n",
        "$$\n",
        "\n",
        "where $\\nabla_\\theta^2 F$ denotes the Hessian matrix, whose $(i, j)$ entry is $\\frac{\\partial^2 F}{\\partial \\theta_i \\partial \\theta_j}$.\n",
        "\n",
        "The left-hand side of this equation is called the **Fisher information** of the family $\\left\\{p(\\cdot | \\theta)\\right\\}_{\\theta \\in \\mathcal{T}}$ at parameter vector $\\theta_0$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "KF-ac0Bk-HmR"
      },
      "source": [
        "#### Proof of claim\n",
        "\n",
        "We have\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "\\left(\\nabla_\\theta^2 \\log p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "\\right]\n",
        "\u0026\\stackrel{\\text{(1)}}{=} \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "  \\left(\\nabla_\\theta^\\top \\frac{\n",
        "    \\nabla_\\theta p(Y | \\theta)\n",
        "  }{\n",
        "    p(Y|\\theta)\n",
        "  }\\right)_{\\theta=\\theta_0}\n",
        "\\right] \\\\\n",
        "\u0026\\stackrel{\\text{(2)}}{=} \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "  \\frac{\n",
        "    \\left(\\nabla^2_\\theta p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\n",
        "  -\n",
        "  \\left(\\frac{\n",
        "    \\left(\\nabla_\\theta\\, p(Y|\\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\\right)\n",
        "  \\left(\\frac{\n",
        "    \\left(\\nabla_\\theta\\, p(Y|\\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\\right)^\\top\n",
        "\\right] \\\\\n",
        "\u0026\\stackrel{\\text{(3)}}{=} \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "  \\frac{\n",
        "    \\left(\\nabla^2_\\theta p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\n",
        "  -\n",
        "  \\text{score}(Y, \\theta_0)\n",
        "  \\,\\text{score}(Y, \\theta_0)^\\top\n",
        "\\right],\n",
        "\\end{align*}\n",
        "\n",
        "where we have used (1) chain rule for differentiation, (2) quotient rule for differentiation, (3) chain rule again, in reverse.\n",
        "\n",
        "To complete the proof, it suffices to show that\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "  \\frac{\n",
        "    \\left(\\nabla^2_\\theta p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\n",
        "\\right]\n",
        "\\stackrel{\\text{?}}{=}\n",
        "0.\n",
        "$$\n",
        "\n",
        "To do that, we pass differentiation under the integral sign twice:\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)}\\left[\n",
        "  \\frac{\n",
        "    \\left(\\nabla^2_\\theta p(Y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(Y|\\theta=\\theta_0)\n",
        "  }\n",
        "\\right]\n",
        "\u0026= \\int_{\\mathcal{Y}}\n",
        "  \\left[\n",
        "  \\frac{\n",
        "    \\left(\\nabla^2_\\theta p(y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  }{\n",
        "    p(y|\\theta=\\theta_0)\n",
        "  }\n",
        "  \\right]\n",
        "  \\, p(y | \\theta=\\theta_0)\\, dy \\\\\n",
        "\u0026= \\int_{\\mathcal{Y}}\n",
        "  \\left(\\nabla^2_\\theta p(y | \\theta)\\right)_{\\theta=\\theta_0}\n",
        "  \\, dy \\\\\n",
        "\u0026= \\left[\n",
        "    \\nabla_\\theta^2 \\left(\n",
        "      \\int_{\\mathcal{Y}} p(y | \\theta) \\, dy\n",
        "    \\right)\n",
        "  \\right]_{\\theta=\\theta_0} \\\\\n",
        "\u0026= \\left[\n",
        "    \\nabla_\\theta^2 \\, 1\n",
        "  \\right]_{\\theta=\\theta_0} \\\\\n",
        "\u0026= 0.\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "kAIJfX7IX_lP"
      },
      "source": [
        "### Lemma about the derivative of the log partition function\n",
        "\n",
        "If $a$, $b$ and $c$ are scalar-valued functions, $c$ twice differentiable, such that the family of distributions $\\left\\{p(\\cdot | \\theta)\\right\\}_{\\theta \\in \\mathcal{T}}$ defined by\n",
        "\n",
        "$$\n",
        "p(y|\\theta) = a(y) \\exp\\left(b(y)\\, \\theta - c(\\theta)\\right)\n",
        "$$\n",
        "\n",
        "satisfies the mild regularity conditions that permit passing differentiation with respect to $\\theta$ under an integral with respect to $y$, then\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ b(Y) \\right]\n",
        "= c'(\\theta_0)\n",
        "$$\n",
        "\n",
        "and\n",
        "\n",
        "$$\n",
        "\\text{Var}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ b(Y) \\right]\n",
        "= c''(\\theta_0).\n",
        "$$\n",
        "\n",
        "(Here $'$ denotes differentiation, so $c'$ and $c''$ are the first and second derivatives of $c$. )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "CYBH-KwpfWhr"
      },
      "source": [
        "#### Proof\n",
        "For this family of distributions, we have $\\text{score}(y, \\theta_0) = b(y) - c'(\\theta_0)$.  The first equation then follows from the fact that $\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ \\text{score}(y, \\theta_0) \\right] = 0$.  Next, we have\n",
        "\n",
        "\\begin{align*}\n",
        "\\text{Var}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ b(Y) \\right]\n",
        "\u0026= \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ \\left(b(Y) - c'(\\theta_0)\\right)^2 \\right] \\\\\n",
        "\u0026= \\text{the one entry of } \\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ \\text{score}(y, \\theta_0) \\text{score}(y, \\theta_0)^\\top \\right] \\\\\n",
        "\u0026= \\text{the one entry of } -\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ \\left(\\nabla_\\theta^2 \\log p(\\cdot | \\theta)\\right)_{\\theta=\\theta_0} \\right] \\\\\n",
        "\u0026= -\\mathbb{E}_{Y \\sim p(\\cdot | \\theta=\\theta_0)} \\left[ -c''(\\theta_0) \\right] \\\\\n",
        "\u0026= c''(\\theta_0).\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "AYpWUvvKcX-e"
      },
      "source": [
        "## Overdispersed Exponential Family\n",
        "\n",
        "A (scalar) **overdispersed exponential family** is a family of distributions whose densities take the form\n",
        "\n",
        "$$\n",
        "p_{\\text{OEF}(m,  T)}(y\\, |\\, \\theta, \\phi) = m(y, \\phi) \\exp\\left(\\frac{\\theta\\, T(y) - A(\\theta)}{\\phi}\\right),\n",
        "$$\n",
        "\n",
        "where $m$ and $T$ are known scalar-valued functions, and $\\theta$ and $\\phi$ are scalar parameters.\n",
        "\n",
        "*\\[Note that $A$ is overdetermined: for any $\\phi_0$, the function $A$ is completely determined by the constraint that\n",
        "$\\int p_{\\text{OEF}(m,  T)}(y\\ |\\ \\theta, \\phi=\\phi_0)\\, dy = 1$\n",
        "for all $\\theta$.  The $A$'s produced by different values of $\\phi_0$ must all be the same, which places a constraint on the functions $m$ and $T$.\\]*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "IgpoijwPf7TV"
      },
      "source": [
        "### Mean and variance of the sufficient statistic\n",
        "\n",
        "Under the same conditions as \"Lemma about the derivative of the log partition function,\" we have\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta, \\phi)} \\left[\n",
        "T(Y)\n",
        "\\right]\n",
        "=\n",
        "A'(\\theta)\n",
        "$$\n",
        "\n",
        "and\n",
        "\n",
        "$$\n",
        "\\text{Var}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta, \\phi)} \\left[\n",
        "T(Y)\n",
        "\\right]\n",
        "=\n",
        "\\phi A''(\\theta).\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "gyf51flphGOK"
      },
      "source": [
        "#### Proof\n",
        "\n",
        "By \"Lemma about the derivative of the log partition function,\" we have\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta, \\phi)} \\left[\n",
        "\\frac{T(Y)}{\\phi}\n",
        "\\right]\n",
        "=\n",
        "\\frac{A'(\\theta)}{\\phi}\n",
        "$$\n",
        "\n",
        "and\n",
        "\n",
        "$$\n",
        "\\text{Var}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta, \\phi)} \\left[\n",
        "\\frac{T(Y)}{\\phi}\n",
        "\\right]\n",
        "=\n",
        "\\frac{A''(\\theta)}{\\phi}.\n",
        "$$\n",
        "\n",
        "The result then follows from the fact that expectation is linear ($\\mathbb{E}[aX] = a\\mathbb{E}[X]$) and variance is degree-2 homogeneous ($\\text{Var}[aX] = a^2 \\,\\text{Var}[X]$)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "mYOnAZv9d4XH"
      },
      "source": [
        "## Generalized Linear Model\n",
        "\n",
        "In a generalized linear model, a predictive distribution for the response variable $Y$ is associated with a vector of observed predictors $x$.  The distribution is a member of an overdispersed exponential family, and the parameter $\\theta$ is replaced by $h(\\eta)$ where $h$ is a known function, $\\eta := x^\\top \\beta$ is the so-called **linear response**, and $\\beta$ is a vector of parameters (regression coefficients) to be learned.  In general the dispersion parameter $\\phi$ could be learned too, but in our setup we will treat $\\phi$ as known.  So our setup is\n",
        "\n",
        "$$\n",
        "Y \\sim p_{\\text{OEF}(m, T)}(\\cdot\\, |\\, \\theta = h(\\eta), \\phi)\n",
        "$$\n",
        "\n",
        "where the model structure is characterized by the distribution $p_{\\text{OEF}(m, T)}$ and the function $h$ which converts linear response to parameters.\n",
        "\n",
        "Traditionally, the mapping from linear response $\\eta$ to mean $\\mu := \\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot\\, |\\, \\theta = h(\\eta), \\phi)}\\left[ Y\\right]$ is denoted\n",
        "\n",
        "$$\n",
        "\\mu = g^{-1}(\\eta).\n",
        "$$\n",
        "\n",
        "This mapping is required to be one-to-one, and its inverse, $g$, is called the **link function** for this GLM.  Typically, one describes a GLM by naming its link function and its family of distributions -- e.g., a \"GLM with Bernoulli distribution and logit link function\" (also known as a logistic regression model).  In order to fully characterize the GLM, the function $h$ must also be specified.  If $h$ is the identity, then $g$ is said to be the **canonical link function**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "t-mrWHH2-wtv"
      },
      "source": [
        "### Claim: Expressing $h'$ in terms of the sufficient statistic\n",
        "\n",
        "Define\n",
        "\n",
        "$$\n",
        "{\\text{Mean}_T}(\\eta)\n",
        ":=\n",
        "\\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(\\eta), \\phi)} \\left[\n",
        "  T(Y)\n",
        "\\right]\n",
        "$$\n",
        "\n",
        "and\n",
        "\n",
        "$$\n",
        "{\\text{Var}_T}(\\eta)\n",
        ":=\n",
        "\\text{Var}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(\\eta), \\phi)} \\left[\n",
        "  T(Y)\n",
        "\\right].\n",
        "$$\n",
        "\n",
        "Then we have\n",
        "\n",
        "$$\n",
        "h'(\\eta) = \\frac{\\phi\\, {\\text{Mean}_T}'(\\eta)}{{\\text{Var}_T}(\\eta)}.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "z36iGKlf_-3F"
      },
      "source": [
        "#### Proof\n",
        "By \"Mean and variance of the sufficient statistic,\" we have\n",
        "\n",
        "$$\n",
        "{\\text{Mean}_T}(\\eta) = A'(h(\\eta)).\n",
        "$$\n",
        "\n",
        "Differentiating with the chain rule, we obtain\n",
        "$$\n",
        "{\\text{Mean}_T}'(\\eta) = A''(h(\\eta))\\, h'(\\eta),\n",
        "$$\n",
        "\n",
        "and by \"Mean and variance of the sufficient statistic,\"\n",
        "\n",
        "$$\n",
        "\\cdots = \\frac{1}{\\phi} {\\text{Var}_T}(\\eta)\\ h'(\\eta).\n",
        "$$\n",
        "\n",
        "The conclusion follows."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "D8LV_QHPx-wV"
      },
      "source": [
        "## Fitting GLM Parameters to Data\n",
        "\n",
        "The properties derived above lend themselves very well to fitting GLM parameters $\\beta$ to a data set.  Quasi-Newton methods such as Fisher scoring rely on the gradient of the log likelihood and the Fisher information, which we now show can be computed especially efficiently for a GLM.\n",
        "\n",
        "Suppose we have observed predictor vectors $x_i$ and associated scalar responses $y_i$.  In matrix form, we'll say we have observed predictors $\\mathbf{x}$ and response $\\mathbf{y}$, where $\\mathbf{x}$ is the matrix whose $i$th row is $x_i^\\top$ and $\\mathbf{y}$ is the vector whose $i$th element is $y_i$.  The log likelihood of parameters $\\beta$ is then\n",
        "\n",
        "$$\n",
        "\\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{y}) = \\sum_{i=1}^{N} \\log p_{\\text{OEF}(m, T)}(y_i\\, |\\, \\theta = h(x_i^\\top \\beta), \\phi).\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "aghNxiO_HFW1"
      },
      "source": [
        "### For a single data sample\n",
        "To simplify the notation, let's first consider the case of a single data point, $N=1$; then we will extend to the general case by additivity.\n",
        "\n",
        "#### Gradient\n",
        "We have\n",
        "\n",
        "\\begin{align*}\n",
        "\\ell(\\beta\\, ;\\, x, y)\n",
        "\u0026= \\log p_{\\text{OEF}(m, T)}(y\\, |\\, \\theta = h(x^\\top \\beta), \\phi) \\\\\n",
        "\u0026= \\log m(y, \\phi) + \\frac{\\theta\\, T(y) - A(\\theta)}{\\phi}, \\quad\\text{where}\\  \\theta = h(x^\\top \\beta).\n",
        "\\end{align*}\n",
        "\n",
        "Hence by the chain rule,\n",
        "\n",
        "$$\n",
        "\\nabla_\\beta \\ell(\\beta\\, ; \\, x, y) = \\frac{T(y) - A'(\\theta)}{\\phi}\\, h'(x^\\top \\beta)\\, x.\n",
        "$$\n",
        "\n",
        "Separately, by \"Mean and variance of the sufficient statistic,\" we have $A'(\\theta) = {\\text{Mean}_T}(x^\\top \\beta)$.  Hence, by \"Claim: Expressing $h'$ in terms of the sufficient statistic,\" we have\n",
        "\n",
        "$$\n",
        "\\cdots =\n",
        "  \\left(T(y) - {\\text{Mean}_T}(x^\\top \\beta)\\right)\n",
        "  \\frac{{\\text{Mean}_T}'(x^\\top \\beta)}{{\\text{Var}_T}(x^\\top \\beta)}\n",
        "  \\,x.\n",
        "$$\n",
        "\n",
        "#### Hessian\n",
        "Differentiating a second time, by the product rule we obtain\n",
        "\n",
        "\\begin{align*}\n",
        "\\nabla_\\beta^2 \\ell(\\beta\\, ;\\, x, y)\n",
        "\u0026=\n",
        "  \\left[\n",
        "    -A''(h(x^\\top \\beta))\\, h'(x^\\top \\beta)\n",
        "  \\right]\n",
        "  h'(x^\\top \\beta)\\, x x^\\top\n",
        "  +\n",
        "  \\left[\n",
        "    T(y) - A'(h(x^\\top \\beta))\n",
        "  \\right]\n",
        "  h''(x^\\top \\beta)\\, xx^\\top\n",
        "  ] \\\\\n",
        "\u0026=\n",
        "  \\left(\n",
        "    -{\\text{Mean}_T}'(x^\\top \\beta)\\, h'(x^\\top \\beta)\n",
        "    +\n",
        "    \\left[T(y) - A'(h(x^\\top \\beta))\\right]\n",
        "  \\right)\\, x x^\\top.\n",
        "\\end{align*}\n",
        "\n",
        "#### Fisher information\n",
        "By \"Mean and variance of the sufficient statistic,\" we have\n",
        "\n",
        "$$\n",
        "\\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x^\\top \\beta), \\phi)} \\left[\n",
        "T(y) - A'(h(x^\\top \\beta))\n",
        "\\right] = 0.\n",
        "$$\n",
        "\n",
        "Hence\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbb{E}_{Y \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x^\\top \\beta), \\phi)} \\left[\n",
        "  \\nabla_\\beta^2 \\ell(\\beta\\, ;\\, x, y)\n",
        "\\right]\n",
        "\u0026=\n",
        "  -{\\text{Mean}_T}'(x^\\top \\beta)\\, h'(x^\\top \\beta) x x^\\top \\\\\n",
        "\u0026=\n",
        "  -\\frac{\\phi\\, {\\text{Mean}_T}'(x^\\top \\beta)^2}{{\\text{Var}_T}(x^\\top \\beta)}\\, x x^\\top.\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "BrA1A583HOng"
      },
      "source": [
        "### For multiple data samples\n",
        "\n",
        "We now extend the $N=1$ case to the general case.  Let $\\boldsymbol{\\eta} := \\mathbf{x} \\beta$ denote the vector whose $i$th coordinate is the linear response from the $i$th data sample.  Let $\\mathbf{T}$ (resp. ${\\textbf{Mean}_T}$, resp. ${\\textbf{Var}_T}$) denote the broadcasted (vectorized) function which applies the scalar-valued function $T$ (resp. ${\\text{Mean}_T}$, resp. ${\\text{Var}_T}$) to each coordinate.  Then we have\n",
        "\n",
        "\\begin{align*}\n",
        "\\nabla_\\beta \\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{y})\n",
        "\u0026= \\sum_{i=1}^{N} \\nabla_\\beta \\ell(\\beta\\, ;\\, x_i, y_i) \\\\\n",
        "\u0026= \\sum_{i=1}^{N}\n",
        "  \\left(T(y) - {\\text{Mean}_T}(x_i^\\top \\beta)\\right)\n",
        "  \\frac{{\\text{Mean}_T}'(x_i^\\top \\beta)}{{\\text{Var}_T}(x_i^\\top \\beta)}\n",
        "  \\, x_i \\\\\n",
        "\u0026=\n",
        "  \\mathbf{x}^\\top\n",
        "  \\,\\text{diag}\\left(\\frac{\n",
        "      {\\textbf{Mean}_T}'(\\mathbf{x} \\beta)\n",
        "    }{\n",
        "      {\\textbf{Var}_T}(\\mathbf{x} \\beta)\n",
        "    }\\right)\n",
        "  \\left(\\mathbf{T}(\\mathbf{y}) - {\\textbf{Mean}_T}(\\mathbf{x} \\beta)\\right) \\\\\n",
        "\\end{align*}\n",
        "\n",
        "and\n",
        "\n",
        "\\begin{align*}\n",
        "\\mathbb{E}_{Y_i \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x_i^\\top \\beta), \\phi)} \\left[\n",
        "  \\nabla_\\beta^2 \\ell(\\beta\\, ;\\, \\mathbf{x}, \\mathbf{Y})\n",
        "\\right]\n",
        "\u0026= \\sum_{i=1}^{N} \\mathbb{E}_{Y_i \\sim p_{\\text{OEF}(m, T)}(\\cdot | \\theta = h(x_i^\\top \\beta), \\phi)} \\left[\n",
        "  \\nabla_\\beta^2 \\ell(\\beta\\, ;\\, x_i, Y_i)\n",
        "\\right] \\\\\n",
        "\u0026= \\sum_{i=1}^{N}\n",
        "  -\\frac{\\phi\\, {\\text{Mean}_T}'(x_i^\\top \\beta)^2}{{\\text{Var}_T}(x_i^\\top \\beta)}\\, x_i x_i^\\top \\\\\n",
        "\u0026=\n",
        "  -\\mathbf{x}^\\top\n",
        "  \\,\\text{diag}\\left(\n",
        "    \\frac{\n",
        "      \\phi\\, {\\textbf{Mean}_T}'(\\mathbf{x} \\beta)^2\n",
        "    }{\n",
        "      {\\textbf{Var}_T}(\\mathbf{x} \\beta)\n",
        "    }\\right)\\,\n",
        "  \\mathbf{x},\n",
        "\\end{align*}\n",
        "\n",
        "where the fractions denote element-wise division."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "jUrOmdt395hZ"
      },
      "source": [
        "## Verifying the Formulas Numerically"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "WVp59IBW-TK6"
      },
      "source": [
        "We now verify the above formula for gradient of the log likelihood numerically using `tf.gradients`, and verify the formula for Fisher information with a Monte Carlo estimate using `tf.hessians`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 0,
      "metadata": {
        "colab": {
          "height": 509
        },
        "colab_type": "code",
        "executionInfo": {
          "elapsed": 127317,
          "status": "ok",
          "timestamp": 1538272038898,
          "user": {
            "displayName": "",
            "photoUrl": "",
            "userId": ""
          },
          "user_tz": 420
        },
        "id": "oM-HDPdPepE-",
        "outputId": "5d311920-7345-4fda-d83d-50104d056eed"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Coordinatewise relative error between naively computed gradients and formula-based gradients (should be zero):\n",
            "[[ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]\n",
            " [ -1.96118673e-16  -3.13789877e-16  -1.96118673e-16]\n",
            " [ -2.08845965e-16  -1.67076772e-16  -2.08845965e-16]]\n",
            "\n",
            "Coordinatewise relative error between average of naively computed Hessian and formula-based FIM (should approach zero as num_trials -\u003e infinity):\n",
            "[[ 0.00072369  0.00072369  0.00072369]\n",
            " [ 0.00072369  0.00072369  0.00072369]\n",
            " [ 0.00072369  0.00072369  0.00072369]]\n",
            "\n"
          ]
        }
      ],
      "source": [
        "def VerifyGradientAndFIM():\n",
        "  model = tfp.glm.BernoulliNormalCDF()\n",
        "  model_matrix = np.array([[1., 5, -2],\n",
        "                           [8, -1, 8]])\n",
        "\n",
        "  def _naive_grad_and_hessian_loss_fn(x, response):\n",
        "    # Computes gradient and Hessian of negative log likelihood using autodiff.\n",
        "    predicted_linear_response = tf.squeeze(\n",
        "        tf.matmul(model_matrix, x[..., tf.newaxis]), axis=[-1])\n",
        "    log_probs = model.log_prob(response, predicted_linear_response)\n",
        "    grad_loss = tf.gradients(-log_probs, [x])[0]\n",
        "    hessian_loss = tf.hessians(-log_probs, [x])[0]\n",
        "    return [grad_loss, hessian_loss]\n",
        "\n",
        "  def _grad_neg_log_likelihood_and_fim_fn(x, response):\n",
        "    # Computes gradient of negative log likelihood and Fisher information matrix\n",
        "    # using the formulas above.\n",
        "    predicted_linear_response = tf.squeeze(\n",
        "        tf.matmul(model_matrix, x[..., tf.newaxis]), axis=[-1])\n",
        "    mean, variance, grad_mean = model(predicted_linear_response)\n",
        "\n",
        "    v = (response - mean) * grad_mean / variance\n",
        "    grad_log_likelihood = tf.squeeze(\n",
        "        tf.matmul(model_matrix, v[..., tf.newaxis], adjoint_a=True), axis=[-1])\n",
        "\n",
        "    w = grad_mean**2 / variance\n",
        "\n",
        "    fisher_info = tf.matmul(\n",
        "        model_matrix,\n",
        "        tf.matmul(tf.diag(w), model_matrix),\n",
        "        adjoint_a=True)\n",
        "    return [-grad_log_likelihood, fisher_info]\n",
        "\n",
        "  # Monte Carlo estimate of E[Hessian(-LogLikelihood)], where the expectation is\n",
        "  # as written in \"Claim (Fisher information)\" above.\n",
        "  num_trials = 20\n",
        "  trial_outputs = []\n",
        "  np.random.seed(10)\n",
        "  model_coefficients_ = np.random.random(size=(model_matrix.shape[1],))\n",
        "  model_coefficients = tf.convert_to_tensor(model_coefficients_)\n",
        "  for _ in range(num_trials):\n",
        "    # Sample from the distribution of `model`\n",
        "    response = np.random.binomial(\n",
        "        1,\n",
        "        scipy.stats.norm().cdf(np.matmul(model_matrix, model_coefficients_))\n",
        "    ).astype(np.float64)\n",
        "    trial_outputs.append(\n",
        "        list(_naive_grad_and_hessian_loss_fn(model_coefficients, response)) +\n",
        "        list(\n",
        "            _grad_neg_log_likelihood_and_fim_fn(model_coefficients, response))\n",
        "    )\n",
        "\n",
        "  naive_grads = tf.stack(\n",
        "      list(naive_grad for [naive_grad, _, _, _] in trial_outputs), axis=0)\n",
        "  fancy_grads = tf.stack(\n",
        "      list(fancy_grad for [_, _, fancy_grad, _] in trial_outputs), axis=0)\n",
        "\n",
        "  average_hess = tf.reduce_mean(tf.stack(\n",
        "      list(hess for [_, hess, _, _] in trial_outputs), axis=0), axis=0)\n",
        "  [_, _, _, fisher_info] = trial_outputs[0]\n",
        "  \n",
        "  with tf.Session() as sess:\n",
        "    naive_grads_, fancy_grads_, average_hess_, fisher_info_ = sess.run(\n",
        "        [naive_grads, fancy_grads, average_hess, fisher_info])\n",
        "  \n",
        "    print(\"Coordinatewise relative error between naively computed gradients and\"\n",
        "          \" formula-based gradients (should be zero):\\n{}\\n\".format(\n",
        "              (naive_grads_ - fancy_grads_) / naive_grads_))\n",
        "    \n",
        "    print(\"Coordinatewise relative error between average of naively computed\"\n",
        "          \" Hessian and formula-based FIM (should approach zero as num_trials\"\n",
        "          \" -\u003e infinity):\\n{}\\n\".format(\n",
        "              (average_hess_ - fisher_info_) / average_hess_))\n",
        "    \n",
        "VerifyGradientAndFIM()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "bAiNubQ-WDHN"
      },
      "source": [
        "# References\n",
        "\n",
        "\u003ca name='1'\u003e\u003c/a\u003e[1]: Guo-Xun Yuan, Chia-Hua Ho and Chih-Jen Lin. An Improved GLMNET for L1-regularized Logistic Regression. _Journal of Machine Learning Research_, 13, 2012.\n",
        "http://www.jmlr.org/papers/volume13/yuan12a/yuan12a.pdf\n",
        "\n",
        "\u003ca name='2'\u003e\u003c/a\u003e[2]: skd. Derivation of Soft Thresholding Operator.  2018.\n",
        "https://math.stackexchange.com/q/511106\n",
        "\n",
        "\u003ca name='3'\u003e\u003c/a\u003e[3]: Wikipedia Contributors. Proximal gradient methods for learning. _Wikipedia, The Free Encyclopedia_, 2018.\n",
        "https://en.wikipedia.org/wiki/Proximal_gradient_methods_for_learning\n",
        "\n",
        "\u003ca name='4'\u003e\u003c/a\u003e[4]: Yao-Liang Yu. The Proximity Operator.\n",
        "https://www.cs.cmu.edu/~suvrit/teach/yaoliang_proximity.pdf"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "Generalized_Linear_Models.ipynb",
      "provenance": [
        {
          "file_id": "1yGEzxU4fWzttyajAPYvbFPwxfZCturhd",
          "timestamp": 1536111336137
        }
      ],
      "version": "0.3.2"
    },
    "kernelspec": {
      "display_name": "Python 2",
      "name": "python2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
